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Abstract 

This work is an analysis of the fluid behavior ina 
pipe subject to a harmonically varying pressure gradient. 
From this basis the acoustic properties are obtained and 
these are compared to experimental data. 

This analysis considers the affect of the inertia, the 
viscosity and the pipe configuration on the fluid motion. 
Galerkin's method of weighted residuals is used to derive a 
Finite Element model from a simplified Navier-Stokes 
momentum equation. A four node isoparametric element is used 
to determine the velocity distribution at two times of the 
harmonic cycle, at the beginning of the cycle (when pressure 
gradient is a maximum), and at 90° into the cycle (when 
pressure gradient is going through zero). Using the mean of 
these velocities the resistivity and the effective density 
are found as functions of frequency. 

For a pipe of uniform cross-sectional area the 
impedance is directly found from the resistivity and the 
effective density via an exact equation. To find the 
impedance for a pipe with an area varying along the axis a 
second Finite Element is developed which requires these same 
parameters. 

Results are presented for the cases of a circular pipe, 
triangular pipe, Square pipe, two parallel walls and a 
rectangular pipe with an aSpect ratio of five. Good 
agreement is found between the predicted and the 


experimental results. 
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Definition of Symbols 


Py Static value of density. 

Cc Speed of sound in tube. 

Pp, Effective density. 

R Resistivity. 

w Circular frequency (rad/sec). 

k | Wave number (w/C). 

u Dynamic viscosity. 

es Hydraulic radius (2area/perimeter). 

K Dimensionless frequency r, (0 w/u)/”, 

U Particle velocity, in X direction. 

U Mean velocity across pipe. 

U Particle acceleration, in X direction. 

tae} Column matrix. 

(eu Rectangular matrix. 

5 Imaginary unity. 

P Maximum acoustic pressure. 
Argument of fluid velocity. 

Rf Reflection factor. | 

|RE | Magnitude of reflection factor. 

Z Specific impedence. 

3 Phase angle. 

aE Shear stress. : 

F Biotis:-[jcorrection: function. 

Y Ratio of specific heats for a gas. 

Ni Interpolation function. 


SEL A: 


| : st 7 ‘ 7 
atoidan! © 026d Sige: «eS 
Neate TO eh aay, sheer’. Ae 


4 ae o oe : " 
oe Gene eer i benge 


an VSI stab. ay ecse Te 


>: : ; 
) Suns - 
+i =” ' hat. Bee | 4 ~< 
| 7 
." | , 7 
4 
, - ow lute 
De ry ‘i cesT (a t oa per LAS PN 
7 ara it owt 
b 
vsiescssns a a 
= » | 
ay AsS AU tee ath ‘3 
a yex 
ie be - : 71190 08 
j 
r » 4a 
, a a | ae) 
1 ' I 
s/ 
: 4 s im ' [en i 
s) 


’ 

i = > a4 © + \ : 
- Bee a A 

, Vie 


: > 
fate > ugae 1 4i¢t 
i a ng ; : J = 
a a? a : 7 Lb pA aus 7 
. ay a | 
a Peri oe ee a a ii he rdeal> *% 
: 


: ate, ~ 
! . f 
i Pe : =4>.: nN : - 7 
Per ih i. 3 the: ‘ 


wee 
= ae ae 

eS 
aie 
» # 7 7 | 
on v2 


Cy 


Domain 
S 
Re(u) 


Im(u) 


Parameters of the field variable, in this work is 
either the velocity at a point or its derivative. 
Entire pipe area. 

Boundary enclosing the pipe area. 

Real part of the complex variable u. 


Imaginary part of the complex variable u. 
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1. Introduction 


The problem to be considered is the analysis of the 
acoustic properties of pipes of arbitrary cross-section in 
which viscous and inertial forces of the fluid are included. 
Previous work in this area has been restricted by the 
analytical methods to the cases of channels and cylindrical 
pipes. 

A review given by Zwikker and Kosten [19] covers and 
builds: on thercontributions of Helmholtz, Kirchhoff. [11]: and 
Rayleigh [15] in the analysis of fluid motion in pipes. The 
first analysis of unsteady fluid motion near a wall was 
developed by Stokes and this analysis was applied to the 
cases of channels and pipes by Helmholtz who considered only 
the viscous and inertial forces. Kirchhoff [11] made a major 
contribution by including viscosity and effects arising from 
the generation of heat and its conduction to and from the 
walls of the tube. Rayleigh [15] built on these works to 
develop a theory of the fluid motion for cases when the 
viscosity dominates and when viscosity 1s very small. 
Rayleigh [15] was the first to suggest that the acoustic 
properties of a wall made of numerous fine channels could 
simulate the properties of a porous sound absorbing wall. 

The seneraimeheractertstics of fluid motion in round 
pipes can best be described through examining the extreme 


cases at low and high frequencies. The most important 
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difference is whether the viscous or the inertial forces 
dominate. 

At low frequencies or when the tube is very small the 
viscous property of the fluid determines its motion. The 
characteristics of the motion are as follows. The mean 
velocity of the fluid is proportional to and in phase with 
the pressure gradient. A boundary condition of viscous flow 
is that the velocity at the pipe wall is zero. At these low 
frequencies in a round pipe the velocity profile is 
approximately a paraboliod of revolution. The temperature of 
the fluid is held constant by the conduction of heat to the 
walls of the pipe, which in most cases has a much larger 
heat capacity than the fluid. (The heat here iene toned: 
refers to that which 1s generated by the compression and 
expansion of the fluid). This results in the speed of sound 
propagation in the pipes being at its lower limit which is 
290 =m/s for arr 

At high frequencies the fluid motion is complicated by 
two factors. The first and most important factor is the 
inertia of the fluid. When the inertial forces dominate over 
the viscous forces the fluid motion lags behind the pressure 
gradient. Under harmonically varying pressure gradients, as 
frequency increases, the motion lag approaches 90°. In 
addition the velocity profile flattens and, for a given 
pressure gradient, the peak velocities achieved decrease 
with increasing frequency. This situation is similar to a 


damper-mass system; if the forcing function is kept 
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constant, the amplitude of vibration decreases when the 
frequency increases. 

A more involved factor influencing the fluid motion at 
high frequencies is the conduction of heat. At high 
frequencies the fluid does not communicate heat fast enough 
to the pipe walls to make reasonable the assumption of 
constant temperature throughout the fluid. This in itself 
has three effects. The first effect is that the speed of 
sound approaches the adiabatic speed of sound which is 343 
m/s for air. Since the walls of the pipe have a decreased 
influence on the temperature there exists a Significant 
temperature difference between the walls and the fluid at 
the axis of the pipe. Thus the viscosity is no longer 
constant. The third effect is the introduction of another 
mechanism for dissipation in the form of thermal processes. 
Even at high frequencies, however, the dissipation of the 
wave energy 1S dominated by viscous forces. 

Biot [1] develops a function in an attempt to describe 
the interaction of the viscous and the inertial forces in 
narrow channels and round pipes. This function is applicable 
to absorptive material but here it is used mainly as a 
comparison to the numerical analysis. 

Biot's function is related to the momentum balance as 


follows dite ao k= Pe Ue: Fx) Riou: eed 
K=o Kz=0 
A more convenient form with which to analyze the 


acoustic properties is in the equation 
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where U is the complex mean velocity (averaged across the 
Pipe at t=0 and t=90° in the cycle), U is its time 
Betivert ve, @reue called ™the elfective density and ’R/is)the 
resistivity. This relation assumes that the fluid velocity 
is constant across the pipe and equal to U. Kuttruff [12] 
assumes that the effective density and resistivity are 
constant but admits that from Cremer [4] this analysis is 
limited to conditions where y js wR aa 

The proceeding analysis 1S not under this constraint. 

It can be reasoned that the parameters of effective 
density and resistivity must change with frequency by simply 
contrasting the fluid motion characteristics at low and high 
frequencies. For example at high frequencies the analogy to 
the damper-mass system shows that the magnitude of the 
velocity decreases with frequency. This implies that the 
resistivity increases with frequency. These properties, 
which are functions of frequency, are also strongly 
dependent on the type of pipe considered. 

The dependence on the type of pipe considered is due to 
the different fluid motion found in different pipes. The 
simplest example is to consider steady state flow in pipes. 
The shear stresses on the walls of non-circular pipes can be 
extremely variable along the perimeter of the pipe. Thus it 
is reasonable to say that the fluid behavior, and thus the 


acoustic properties, in non-circular pipes at various 
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frequencies are significantly different from those of the 
cases already studied, i.e. the channel and cylinder. 

The aim of this thesis is to develop an analysis of 
fluid motion which is applicable to pipes of any arbitrary 
cross-section. 

Biot's [1] function for the channels and the 
cylindrical pipe is used to confirm the results obtained. 
Biot [1] states that these two cases represent extreme 
cases. So the function for any other case should lie between 
these two functions. It is found that the cases of 
non-circular pipes do not behave quite as predictably as 
this reasoning suggests. 

More importantly, the intent of this analysis is to 
find how the resistivity and effective density vary with 
frequency for a number of shapes of pipe. With these 
relations it is possible to use the work of Zwikker and 
Kosten [19] and that of Craggs [2] to deduce first the input 
impedance, then directly the absorption, of a matrix of 
finite pipes. This done it is possible to confirm the 
analysis by experiment. This was done along with the above 
analysis. 

The presentation of the analysis proceeds as follows. 
In Chapter 2 the method of analysing the problem is covered. 
This begins with the assumptions which lead to the governing 
equations. The method of solution to the governing 


equations, the finite element method, is then developed. 
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Chapter 3 presents the solution of the governing 
equations as applied to triangular pipes and square pipes. 
Comparison is made of the results of the analysis to Biot's 
bth functions: 

Chapter 4 shows the types of meshes used and the 
convergence of the bulk properties due to changing the type 
of mesh. 

The resulting shear stresses on the walls of the 
triangular pipe and square pipe at various frequencies are 
noted in Chapter 5. 

The relationship between the fluid motion and the 
acouStic properties of resistivity and effective density are 
derived in Chapter 6. These properties are presented, as 
predicted by the numerical model, for the cylindrical, 
triangular, square, and rectangular pipes as well as for the 
slit. Also included are two ways of finding the impedance 
for a finite pipe with one end closed. The first way is from 
the exact equation for a pipe which is uniform along its 
length. The second is a finite element method which can 
model a pipe varying in area along its axis. 

The impedance tube experiment, the predicted absorption 
of the numerical model, and the most critical assumptions of 
the analysis are di aeuseed in Chapter 7. 

A summary of the work and the conclusions drawn from it 


are contained in Chapter 8. 
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2. Calculating Velocity Distributions due to Oscillating 


Pressure Gradients 


2.1 Deriving Velocity Profiles in Small pipes 

The first area to be examined is the fluid velocity 
distribution across a pipe with the fluid being subjected to 
an oscillatory pressure gradient. The importance of this is 
to provide the mean velocity as a function of time. This 
information directly leads to two important bulk properties: 


. These 


the resistivity R, and effective fluid density, Po 


two properties depend on the cross-sectional geometry of the 
pipe and frequency of the oscillatory pressure gradient. 

The acoustic impedance can be determined knowing only 
the resistivity, effective density, and the length between 
the open end where the wave enters and the rigid terminus. 
In practice the acoustic impedance and attenuation are 
required for an assembly or panel of parallel pipes of 
Similar cross-sectional geometry which may be used for 
absorbing sound. Each pipe in the assembly affects the 
portion of the sound wave which strikes it independently and 
identically to the surrounding pipes, therefore the 
reflected wave from the single pipe will have the same 
Strength and phase as the reflected wave from the assembly 
of pipes. 

Pipes of circular cross-section are modelled to compare 


with the previously known analytical results. Pipes of 
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triangular and square cross-section are modelled to compare 
with the experimental results. Small pipes of triangular and 
square cross-section are quite practical due to their 


perfect packing density and ease of manufacture. 


2.2 Assumptions 

To find the velocity distribution certain assumptions 
can be made. 

1) The fluid is Newtonian. 

2) The walls of the pipe are rigid and stationary. 

3) The pressure gradient is strictly axial. 

4) The inlet effects are negligible. 

5) The fluid temperature is constant. 

6) The fluid density is constant. 

The inlet effects are neglected because for significant 
sound absorption the pipes must be fairly thin and long to 
provide, in essence, a large area to act on the fluid. For 
steady flow, the section of pipe taken for the flow to go 
from the initially flat velocity profile to its steady state 
parabolic velocity profile (for laminar flow) is 
proportional to the Reynolds number [16]. In oscillating 
flow the inlet effects decrease for the following reasons. 
First of all the Reynolds number (based on instantaneous 
velocity) decreases with increasing frequency, and second 
the velocity distribution changes from a parabolic toa 


flattened profile thus localizing the shear to the walls and 
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thus the pseudo-steady-state condition is more easily 
established. The effect of the inlet on pulsatile fluid 
motion is covered briefly by McDonald [14] . 

A more Significant assumption is that there is no 
temperature gradient in the fluid. This allows the use of a 
constant viscosity across the pipe at all times but it means 
the absorption of the model will be low because the thermal 
process being neglected is a mechanism of dissipation. This 
assumption means that the isothermal speed of sound C=290 
m/s should be used. 

It is also assumed that the density variation is 
negligible. Although this analysis is to be applied to 
peter neste characteristics of sound waves in pipes the 
fluid motion can be found first by assuming the fluid is 
incompressible. This is reasonable since even a sound. 
pressure level of 100 dB creates a density variation which 
is five orders of magnitude smaller than the average 


density. 


2.3 Governing Equations 
With the above assumptions the Navier-Stokes momentum 


equation [16] is 
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For fluid motion ina stationary pipe the boundary 
condition is U=0 at the pipe wall. The Cartesian coordinate 


system is used with the X direction being along the axis of 
h i © 
the pipe y 


a 


The second governing equation is the wave equation, as 
derived in Appendix A. This equation provides the form of 


the pressure variation needed. In one dimension the equation 
is 


The resulting pressure is 


Bey = Prexpt yp Git -7kx)} 23 


where P is the peak pressure, w is the circular frequency, k 
is the wave number k= w /C, and C is the sonic speed. 


The velocity must assume the same form. 


“Cx ,y,z,t) % Yey,z) aN ANOS oe 


The phase difference between the velocity and the pressure 


is expressed by the complex quantity U. The velocity output 
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of the computer program will be in terms of Ur and Ui sawhicn 


“a 


are the real and imaginary parts of U respectively 
U=U_ + 4U, 2.5 
1g 1 


Looking at one position of the pipe, say X=0, and 
replacing the exponential notation with its equivalent gives 
us 

Uton vice anor UL cos (amt) -7SsinGao. t) ) 
eds) 

Our input quantity 1S pressure so P is set to be real 
only 


Pew ipicos (at) Dy 


Since only the real part of these quantities has 


physical significance 


UConvnzetue=-ULacos( aati Ul san iat) 
208 
Ur is the velocity distribution at t=0 and the negative of 
Ui is the velocity distribution at w t=90°. 
| Substituting this harmonic velocity into the momentum 
equation yields 
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2.4 Galerkin' s Method and the Finite leer Ps Method 

Analytical solutions to the above equation have been 
found in the case of a round pipe as well as for two 
infinite parallel plates, referred to hereafter as a slit. 
To solve for the velocity profile, however, in a pipe with 
an arbitrary cross-section, numerical techniques must be 
employed. 

The numerical technique chosen is Galerkin's method ' 
which is a specific form of the weighted residuals method. 
Huebner [10], who provides some excellent examples of the 
use of Galerkin's method in finite elements, offers the 
following which can be taken as a synopsis of the method. 

"Applying the method of weighted residuals involves 
basically two steps. The first is to assume the 
general functional behavior of the dependent field 
variable in some way so as to approximately satisfy 
the given differential equation and boundary 
conditions. Substitution of this approximation into 
the original differential equation and boundary 
conditions then results in some error called a 
residual. This residual is required to vanish in 
some average sense over the entire solution domain". 
"The second step is to solve the equation (or 
equations) resulting from the first step and thereby 
specialize the general functional form to a 


‘Por further discussion see Refs. 6 and 8. The original work 
by Galerkin is not readily accessible but brief accounts are 
contained in various Russian publications eg. Refs. 5 and 9. 


a) 

az 
a Py 
sa rf j 


eer) sae | ae i see re 


gs #A 1 Gisa dieu 2 ee + eden an 


os a ; a } ee © a 1 V4 


—- : 
os Pu ; 
: ; ' i weiss 
; ‘i - , &é 
wf a ry rs fi ANe : ee OO a) weer : 7 
7 ’ 7 } i E : 
es es s targa ee ee 
7 : 70? re Drs ve i ; 
Aon a ae i 
a] 4 2 i : A a’ 
Pathe til i9 BaF Gas £e 


i 
} a _ 
7 . 


4 O00 ES fee} aes rloeks 


. e534 in Say ea io" fs, Cun am 
Lh ‘i me PIc.. ae ty: atte a 


ey 
aly, ied Se SADIE SS 


r4nlia. beat 1 tdi 

jyiO) aes aade79 oe sog Ji naie 
ob 4 bbe Wal RBA (Na tise wed 

L2e Ours Sip ee rN pabtt i Acw 
wes depat Bh La i) eT | Saya 
a. \ fire cedd (yese ohNGs “ORNS bore 
= ae: pulae 6: gs. | uu'prennebed ese" 


a sitet SOG eet a1, ale a a 


pera DBs £6 78r GS mie We. « were a 


@ &= +) tae ai OG em ase 


Pod 
Pete MNS sictheit lh hoe a 
eee son” )'p: 

Mak Weave Sonat oi 


Le 


particular function, which then becomes the 
approximate solution sought". 

Let us expand this with the simplified momentum 
equation (2.9) in minds The picid variable is the velocity 
U. The domain of interest is the YZ plane bounded by the 
Pipe walls, S. 

First U is approximated by u, where either the 
functional behavior of & is specified in terms of unknown 
parameters, or the functional dependence on all but one of 
the independent variables is specified while the functional 
dependence on the remaining independent variable is left 


unspecified. Then the dependent variable is approximated by 


where Ni are the assumed functions and Ci are either the 
unknown parameters or unknown functions of one of the 
independent variables. The m is the number of unknown 
Parameters with which U is approximated. The m functions Ni 
are generally chosen to satisfy the boundary conditions. 
When the approximation U is substituted into the 
momentum equation (2.9) there will likely be some error 


called the residual, e« . 
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The focus of the method of weighted residuals is to 
determine the m unknowns Ci in such a way that the residual 
over the entire domain is small. 

Therefore we choose m linearly independent weighting 
functions Wi and specify that the weighted average of the 


error vanish over the domain. That is if 


j e W, dA = 0 2.12 


Domain 


then the residual is small in some sense. Note that 
minimizing the error of the differential equation does not 
necesSarily mean that the error of the solution has been 
minimized by the same weighting functions. 

At this point a broad choice of weighting functions are 
available. For this problem, however, Galerkin's criterion 
determines the weighting functions. Galerkin's criterion is 
that the weighting functions are the same as the 
approximating functions used to represent U, that is Wi=Ni 


LOVE Aalst 


2.5 Discretization 

So far the entire solution domain has been considered. 
This requires that the approximating functions identically 
Satisfy the boundary conditions. This can mean that the 
functions can be quite involved if the shape of the boundary 


is irregular. Since the differential equation holds at any 
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point in the domain, it also holds for any collection of 
points defining an arbitrary subdomain or element of the 
whole domain. Focusing on this element generally allows a 
Simplification of the approximating function. The shape of 
the domain is then approximated by an assembly of such 
elements. 

Now, looking at a typical element, the functions Ni are 
the interpolation functions defined over the element, and 
the undetermined parameters Ci are the values of the field 
variables or their derivatives at the nodes. For each 


element, then, 


f 28 + jo wu ©? Sts Ui ena ak sD 2.13 
p 62? 
where 
1g 
Cee Nice pc) 2.14 
i=1 | 


where r is the number of unkown parameters of each 

element i.e. the product of the number of nodes per element 
and the number of unknowns per node. The unknowns are 
velocity if the lowest order interpolation function is used 
or the unknowns may also include the derivatives of the 
velocity if a higher order function is used. Henceforth U 
will represent the approximate velocity. 

The only problem remaining before assembling the 
elements and solving the equations is that the interpolation 
function must ensure continuity of the function between 


elements. The interpolation functions must be chosen such 
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that, at the element boundaries, continuity of the field 
variable, as well as continuity of its partial derivatives 
is Satisfied (up to one order less than the highest-order 
derivative appearing in the expression to be integrated). 
The higher the order of continuity required, the fewer 
Suitable interpolation functions available. 

Applying integration by parts to the integral 
expression, equation (2.9), reduces the order of the 
derivatives and the required order of the approximating 
function. Integration by parts also conveniently introduces 
the natural boundary conditions. 

Appendix B shows the detailed derivation of the matrix 
equations. For this work the simplest quadrilateral element 
was used. It has only one degree of freedom at each node and 
thus can only assume a linear variation between nodes. 
Figure 2.1 shows a typical element and the interpolation 


function in the transformed ZY plane. 
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Figure 2.1 Coordinate Systems 
The YZ plane is the physical plane. The}§ plane is the 


natural coordinate system. This transformation allows us 
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to perform the necessary integrations using Gauss quadrature 
in an orthogonal plane and yet allowing the element to 
assume the shape of practically any quadrilateral in the ZY 
plane. The transformation from one plane to the other is 
performed by the same function as the interpolation 
Function, 

Z = +x {Z, (1+) (1-n)+Z, (1-6) (14n) + 


2, (146) (1-n) + Z, (146) (14n)} 2.16 


Y= 2x {¥,(1-€)(1-n) + Z,(1-€)(14n) + 
Yq (1+€) (1-n) + Y, (+8) (14n) yapenly! 
Elements calculated by this procedure are termed 


isoparametric. 


2.6 Matrix Equations 
After applying Galerkin's method to the differential 
equation, equation (2.13) must be solved for each element 


domain. In matrix notation 


@yi) dkseste} 2.18 
aX. 
D 
jou fNi U dA = jpjw [Me] {U} 2.19 
D 
u f Niv? u dA = 4u [Del {Ul 2.20 
D 


Where {P} is an array containing the force of pressure 
assigned to each node, that is, the larger the element the 


larger the values of force at each node. {U} is the array of 
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the velocities at each node in the element, [Me] is called 
the mass matrix, it represents the mass of the element, and 
[De] is the stiffness matrix. So the equation can be written 


for each element as 


{P}= (u[De] + j pW [Me]) {U} 22k 


These elemental equations are assembled into global matrices 
producing simultaneous equations to be solved. 

If a certain procedure is applied to the global 
numbering of the mesh an important process can be applied to 
the global matrices. The process is partitioning. The 
required system is numbering the nodes on the interior of 
the pipe first leaving the boundary nodes to be numbered 
last. The resulting global matrix will have all of the 
information relevent to the interior nodes in the top left 
corner and all of the information relevent to the boundary 
nodes in the bottom left corner. This will provide a direct 
means to finding the shear forces on the pipe walls after 
solving for the velocities at the interior nodes. 

Define: 
DOF = number of degrees of freedom at each node allowed by 
interpolation functions. 


NK DOF times the number of nodes on pipe wall. 


NI = DOF times the number of nodes inside pipe i.e. where u 
is not constrained by the boundary conditions. 


NT = NK + NI 
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Array size Array Comments 
NIx1 {U} = matrix of velocities at interior nodes. 
NKx 1 {Uo} = matrix of velocities at nodes on pipe 


Wall, (filled) with: zeros) . 
NIx1 {P} = force due to the pressure at each 


internal node 


NKx 1 {F} = force due to shear at each node on 
pipe wall. 

NTXNT [Ma] = assembly of all mass matrices [Me]. 

NIxNI [M,,] = subset of [Ma] , represents mass of 


fluid around each interior node 
NKXNI [M._.37 = subset of [Ma] , represents mass of 


fluid around each node on pipe wall. 


NTXNT [Da] = assembly of all stiffness matrices 
[De]. 
NIXNI [D,,] = subset of [Da] | represents. stiffness 


of fluid at each interior pipe node. 


NKxNI [De = subset of [Daly vstifinesssot sfluid at 


Car 
nodes on pipe walls 


The partitioning appears as follows 


[Mal = [Da} = Boe. 
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The result of assembling all the elemental equations is 


{F} Oe eee ne bce je f 
Only viUl sand {F}-are’ unknown. fU} lis selved, from 


{u}- wb] +see(e by iP} jiepy) 


It 1S interesting to see that the shear forces on the 
nodes of the walls can now be obtained by a single matrix 


multiplication 


{F}= @lp,} + ipnwled) fu} DES 


This last step iS made easy pecetce of the aforementioned 
numbering scheme. This could be made without this scheme but 
only if intricate matrix manipulations are made to put the 
resulting matrices into the form already shown above. 

Some computing facilities can not solve systems which 
consist of complex variables, as these do. Therefore it may 
be desirable to express the above in terms of real matrices 
only. The complex equation to find {U} can be written in 


terms of real numbers only as 
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where subscript r denotes the real part and subscript i the 
imaginary part of the quantity at each node. 


Similarly treating the shear force equation results in 


Fr »DeF] — fw [Mce] Ur 
= eon 


Fi pwlMce] kK (pet | Ui 

Thus the stage has been reached where the problem has 
been reduced to determining an appropriate mesh of elements 
and solving sets of linear equations. 

Determining the best element mesh,it will be shown in 
the next chapter, is done interactively with the velocity 
aistribution:. 

The assumptions made and numerical technique employed 
provide a basis on which to obtain the fluid motion and 
ultimately the acoustic characteristics of a pipe. Up to 
this point it is difficult to see what limitations might 
occur under these assumptions and the numerical analysis. On 
examination of the results it can be seen that the analysis 
is limited to certain frequencies. By choosing an 
appropriate finite element mesh, however, the range is 
considerably extended. This will be further discussed in 
chapter 4. 

This leads to the initial results, the velocity 


profiles. 
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3. Resulting Velocity Distributions 


3.1 Interaction of Inertia and Viscosity 

Now consider the properties of the velocity associated 
with an oscillatory pressure gradient. A physical 
description will be given followed by the Finite Element 
results. First of all note a parameter commonly related to 


pulsatile flow problems, K, defined as ' 


Pow ie 
Cri hee 31 


K 1S similar to Reynold's number in that it compares the 
inertial to the viscous forces. Commonly K is called either 
the acoustic Reynold's number or the dimensionless 
frequency. 

The first feature of the fluid motion to be noted is 
that there is a phase lag between the pressure gradient and 
the fluid movement. The fluid near the wall is almost in 
phase with the pressure gradient. This is because the 
viscous losses at the wall prevent this layer from gaining 
much momentum. The less momentum the greater the ease of the 
fluid to follow the pressure gradient reversals. The fluid 
near the axis of the pipe, which is least affected by the 


walls, achieves the highest momentum, and lags the most 


‘ r is some typical dimension. For a slit r is half the 
width for all other cases it is taken as the wetted 
perimeter radius or hydraulic radius defined as twice the 
area divided by the perimeter. 
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behind the pressure reversals. At high values of K there is 
less time in the cycle for the wall shear to translate to 
the centrally located fluid and the velocity profile of the 
interior fluid becomes flattened. The fluid in the central 
part of the tube is virtually unsheared and the very high 
velocity gradients are found near the wall. Thus the fluid 
behaves like a solid core surrounded by a thin viscous fluid 
(See Figure 3.1). According to McDonald [14] at a relatively 
high frequency of K=20 the region of high velocity gradient, 
where the viscosity is of any importance, is virtually 
confined to the outer 5 per cent of the tube radius (See 


Pagure 322): 


3.2 Velocities in Triangular and Square Pipes 

The velocity distributions for the triangular and 
Square pipes are shown for various frequencies in Figures 
Seo tO F1lgures 73 ..i'. 

It is easy to see the Similarity of these flows to 
those of the round pipe, as just described. As the frequency 
increases the velocity distribution in the middle of the 
pipe flattens and becomes more out of phase with the 
pressure gradient. Correspondingly the peaks in the corners 
of the pipe at t=0 become accentuated. The three dimensional 
figures do not have the same velocity scale. To show the 
relative changes of the velocities as the frequency 


increases see Figure 3.12 and Figure 3.13 . Shown here are 
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the velocity profiles for a round pipe with the frequency 
changing but with the pressure gradient held constant. 
Figure 3.13 is the velocity profile at time t=0 in the 
cycle. Figure 3.14 is the velocity profile at time wt=90°, 
note that the scale of both figures is the same. 

It might be noticed that the velocity profiles for the 
triangular and square pipes contain some irregulaties. The 
major cause of these irregularities is the process used for 
plotting. The input data for the plotting program must be 
from a mesh of uniform spacing throughout. This means that 
the output of the F.E programs must be used to interpolate 
values at each grid point for the uniform mesh. These 
interpolations were handled by a system library program 
which unfortunately did not use the same interpolation 
functions as those of the elements. At high values of K the 
irregularities are especially prominant near the edges of 
the pipe. The points that suddenly stick out are points 
where the nodes of the original Finite Element grid points 
are. The velocities in the corners of the pipes are better 
represented because the original grid contains a higher 
concentration of nodes in the corners. Some other 
irregularities originate in the plotting program itself. 
This is seen from the erratic lines drawn along the edge of 
the square pipe even though the input function is zero along 


the edge. 
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Figure 3.1] 


Velocity Profiles in circular 
times in the cycle. 
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Figure 3.2 Velocity profiles of a cylindrical pipe created 
by a pressure gradient of a sinusoidal wave form. 
Profiles are given for half the tube at intervals 
Of -oOee from Osco. L300 
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Figure 3.5 Contours of velocity in triangular pipe. 
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Figure 3.8 Contours of velocity in triangular pipe. 
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Figure 3.11 Contours of velocity in triangular pipe. 
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Figure 3.14 Contours of velocity in triangular pipe. 
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Figure 3.17 Contours of velocity in triangular pipe. 
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Figure 3.18 Velocity profile in triangular pipe. 
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Figure 3.20 Contours of velocity in triangular pipe. 
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Figure 3.39 Circular pipe velocity profile at t=0 for 
various K. Pressure gradient held constant. 


63 


—— a +. fn ~ : 
: = — —- ——— —— 2 


4-0 5.0 - 


yo? O~? Se. o Le eg ge 
ah seeno flu; root he ll ut 
7 7 


4 


0.30 


—ODPN® 


ao 


0.25 


K 
K 
K 
K 
K 
K 


0.16 0.20 


0.10 


~~ 

° 

Ss 

. Ts) 

u = 

key 

py 

4 

ar 

(é) 

as 
o 

cab) 

SA < 
oO 
Ww 
fej 
fam } 
| 


0.0 0.1 O.2 O.8 O.4 O56 O6 O.7 O8 O98 1.0 


Figure 3.40 Circular pipe velocity profile at wt=90° 
for various K. Pressure gradient held constant. 
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4, Convergence and Influence of Mesh type 


4.1 Biot's Function 

The mean values of the real and imaginary velocities 
will be used to arrive at the acoustic properties of the 
Pipe. In order to test the mean values predicted by the 
numerical model some analytical work relating similar mean 
properties to frequency is introduced. It has only been 
possible to apply this work to propagation in a circular 
Pipe and between two parallel plates - heretofore assumed to 
produce functions which would represent bounds on the 
possible functions. The principle results of this section 
are the Biot functions and these are derived primarily to 
verify the results of the finite element program. This 
verification has also provided insight into optimizing the 
convergence of the bulk properties through changing the 
finite element mesh. 


Biot [1] defines a term F(K) (compare to equation 1.1), 


2b 
F epi oned 


(K) Cry 


Gila 


4H is the dynamic viscosity, t is the shear stress;Uis the 
average velocity over a cross-section, (F,t, and uvarevall 
complex, and -fusccions.of K). For a slit b,is half the 


distance between the parallel walls. For a circular pipe, b, 
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is the radius, for the other pipes it is the hydraulic 
radius. 

Biot [1] describes this function F as a measure of the 
deviation from steady flow. Therefore, the constant c is 
introduced. Its purpose is to make the function F equal 
unity at K=0 . The analytical solutions show that for a 
Circular pipe C=8 and for the parallel walls C=6. The 
numerical model solutions find that for a triangular pipe 
C=6.68 for a Square pipe C=7.14 and for a rectangular pipe 
with an aspect ratio of five C=9.59 

The real part of the function F versus K and the 
imaginary part versus K are shown respectively in Figures 
4.1 and 4.2. These graphs compare the curves for the case of 
the parallel walls, the circular, triangular and square 
pipes as well as for a rectangular pipe with an aspect ratio 
eir5s 

AccordingrtosBiot [4)JAfor) thescincular) pipe; the 
functions at high K asymptotically become parallel to 


straight lines 
EieeKe/ av (32) =90al774K: 4.2 
For the slit the curves asymptotically become parallel to 


Eu=yihe/ 1 ¥( 18) =0«236eK. 4.3 
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Biot [1] mentions that the slit and the circular pipe 
should represent the extreme aspect ratios thus the Biot 
curves for all of the other configurations should be between 
these two curves. The results presented here do not appear 
to confirm this. The curves for the rectangle lie just 
Bueside the curve of the circle. This discrepancy may be in 
the definition of a typical dimension to calculate the 
acoustic Reynolds number. At high aspect ratios there is a 
problem using the hydraulic radius as the typical dimension. 
i: As an example of this problem consider the slit. From the 
hydraulic radius equation, the radius is equal to the width 
of the slit whereas for Biot's work half the width is used 
in the calculation of K. This problem has been left 
unresolved. This does not detract from the validity of the 
computer program results since the equations are solved in 
terms of the actual frequency and actual pipe dimensions. 
The only area this influences is the principle of similar 
solutions for the high aspect ratio configurations. 

If the aspect ratio is indeed a major factor 
determining the correction function it would be reasonable 
to expect a triangular pipe and especially a square pipe (as 
its aspect ratio is clearly unity) to have curves closer to 
that of the circular pipe as opposed to that of the slit. 
This conclusion is contradicted by the results in Figure 4.1 
and Figure 4.15. This result is unexpected and can only 


suggest that there is another factor to be considered. 
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4.2 Mesh types 

By comparing the Biot functions resulting from 
different element meshes information can be gleaned to help 
find the characteristics of a mesh which would minimize 
error. 

First the easiest case to model is used to demonstrate 
some of the pap kes used in choosing appropriate meshes. 
This case is the circular pipe. Since the circle is 
axisymmetric a small sector of a circle can be used to model 
the entire pipe. See Figure 4.3 

If an arc of 20° is chosen and 4 quadrilateral elements 
are used to define the sector the result is the same as a 
circle modelled by 72 elements. 

As mentioned, the elements chosen assume linear 
variations and have only one unknown per node. If higher 
order elements are used the use of symmetry iS more involved 
because the derivatives across unconstrained boundaries must 
be specified to be zero. A constrained boundary is wherever 
the velocity is set to zero to satisfy the boundary 
condition. 

The grids shown are selected because, of the grids 
tried, they give the best results for the number of 
unconstrained nodes used. The CPU time required is 
approximately proportional to the square of the number of 
unconstrained nodes. As an example of the time required 
consider the triangular grid (see Figure 4.4) used to model 


the square pipe of Figure 4.6. The mesh uses 74 
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unconstrained nodes. To solve for the velocity distribution 
for a single frequency takes 3.31 CPU seconds on the Amdahl 
470/V8. Most of the CPU is used to invert the matrix of 
equation 2.24. This grid gives the same result as if the 
Square was modelled by 497 nodes. Using the grid of Figure 
4.7 for a rectangular grid finding a single velocity profile 
takes 7.02 CPU sec. 

The CPU time could be reduced by capitalizing on the 
fact that the matrix to be inverted is symmetric. In 
addition certain global numbering schemes allow one to 
minimize the space required by banding the matrix but the 
scheme desired to satisfy the partitioning is not compatible 
with banding. 

Figures 4.8 and 4.9 show that the bulk fluid properties 
are easy to model at low frequencies but a finer finite 
element mesh is required to model the properties at high 
frequencies. It was found that accuracy of a particular 
number of elements could be improved simply by concentrating 
the nodes near the edge of the pipe. More nodes are required 
to model the large changes of the velocity gradients near 
the wall. Large changes are difficult to model in any 
numerical technique. 

To concentrate elements in regions with high gradients, 
the uniform grid in the Z direction is transformed by the 
function 

Znew=Ze exp(1-Z) 4.4 


and the Y values are changed to maintain the straight sides 
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Giut he Sector. 

This new mesh is a tremendous improvement as can be 
seen by comparing Figure 4.8 and Figure 4.9 to Figure 4.10 
and Figure 4.11. For the real part of the Biot function the 
analytical solution is almost indistinguishable from the 
finite element solution with only 4 elements. For the 
imaginary part this change gives each grid the accuracy of a 
uniform grid of roughly 4 times the number of elements of a 
UGVeOrm guid. 

As mentioned the number of elements that effectively 
model a circular pipe depend on the angle of the sector and 
the number of elements in the arc. From this it might be 
“naee nase Surmised that the accuracy can be doubled by 
halving the angle of the sector. If the arc is less than 
about 20° there is, in fact, very little gain in decreasing 
the angle of the arc. The spacing with respect to the pipe 
wall has a much greater significance. By simply looking at 
the velocity profiles some feeling can be acquired as to the 
best places for finer gridding i.e. at the largest changes 
of the velocity gradients. The grids modelling the 
triangular and square pipes contained a higher density of 
elements near the edges and corners of the pipe. 

Figure 4.3 shows the type of grids used to model the 
circle. Velocities at the nodes on the circumference of the 
Girelesare,set touzero, 

Note that for the arc the element nearest the center is 


quite distorted from a square. Large distortions may cause 
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errors in the transformation from the orthogonal plane (in 
natural coordinates) to the global YZ plane. Care must 
be taken when generating meshes that unique mapping from one 
plane to the other takes place. The distortion is acceptable 
if the sign of the Jacobian used in transformations does not 
change over the domain. The distorted element mentioned 
above, by this criterion, is acceptable. Problems may occur 
if one of the interior angles of the quadrilateral is 
greater than 180°. 

Figure 4.4 shows one of the best grids used to model a 
triangular pipe. Its hydraulic radius is one. The horizontal 
axis is Z, the vertical axis is Y. The triangle is modelled 
by a mesh covering only half of the entire domain. The part 
used by the program is shown by Figure 4.4(a). The other 
nodes are just a reflection over the Z axis. The velocities 
are constrained at the nodes on the line Y = Z and Z = 
2.414. 

Figure 4.6 shows how the same mesh is used to model a 
Square. The mesh is scaled down so that the hydraulic radius 
of the square is unity. In this case only the nodes on the 
line Z = 1 are constrained. 

Figure 4.7 shows the effective grid to model the 
rectangle. Appropriate scaling factors in the Y and Z 
directions transformed this mesh into a rectangle with the 
desired aspect ratio and hydraulic radius. This mesh has 121 
nodes and 100 elements. The basic grid is the portion of 


Figure 4.7 of where O<Y<1, 0<ZS1, the constrained portion of 
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the boundary is at Y=1 and Z=1. 

Chapter 6 contains another example and discussion of 
the affect the grid type has on two other bulk properties, 
resistivity and effective density. 

By knowing the type of velocities say from a crude mesh 
a new mesh can be devised with nodes concentrated in the 
appropriate areas, thus minimizing the error in the bulk 
quantities. This is most important because the bulk 
Quantities are needed to derive the acoustic properties of 


the pipe as will be discussed in chapter 6. 


o qolaauottk fhe Sigma 


soicyeenae whee vwdso fp: 


»| ¢ we > yee teh 


a A ; 6128 2 ‘i 


Mb Tm ess Siimandea quds cations tak a 
id #09 eadegud, om (See Bd mh te oe an 
Ps 23 SOA MPO ce: 8a I ‘wh od bgbeatt oa 8 4 'yiga 


0) Seeds ay am bah Lass, + ne wal 


Biot's Function, Re(F) 


166. 9200 9206 3.059355 450 (425 *670 


1.0 


0.0 0.5 


O——-#] SQUARE 
@O@——® CIRCLE 
4———A TRI. 
——x SLIT 
+, (RECS 


Figure 4.1 Real (F) vs. K 
For slit, circle, triangle,square and 
rectangle with aspect ratio of 5. 
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Figure 4.2 Imaginary (F) vs. K. 
For slit, circle, triangle, square, 
rectangle with aspect ratio of 5. 
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(a) For uniformly spaced elements. 
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(b) Improved mesh with 4 elements. 
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(d) Improved mesh with 8 elements. 


Figure 4.3 Sample meshes to model circular pipe. 
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Figure 4.4. (a) Basic 79 node grid to model triangle 
and square. 


(b) Effective grid of triangle using (a) 
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Figure 4.7 Effective grid for square or rectangle 
using 121 node portion in 0<Y<1, 0<Z<1. 
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Figure 4.8 Real (F) vs. K. 
For circle using various uniform grids in arc. 
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Figure 4.9 Imaginary (F) vs. K. 


For circle using various uniform grids in arc. 
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Figure 4.10 Real (F) vs. K. 
For circle using improved grids. 
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Figure 4.11 Imaginary (F) vs. K. 
For circle using improved grids. 
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5. Shear Stresses 
In the derivation of the matrix equations in chapter 2 
it was shown how the forces the fluid exerts on the pipe 
walls could be found. The conventional way of studying the 
forces or stresses would be to use Newton's law' of 


viscosity as given by 


where n is perpendicular to the wall. It is difficult to use 
this equation in this application unless the element meshes 
are chosen such that near the walls there are lines of nodes 
perpendicular to the walls. Even with this type of grid, 
which is inconvenient to generate, errors outa be 
introduced in the process of taking the derivatives of 
velocity. So it is clearly better to use the shear force 
matrix equations as derived in chapter 2. The results that 
this equation produces are the forces in the region of each 
node required to maintain equilibrium in the system. 

A weighting system must be applied to the force ata 
particular node to determine the average shear stress around 
that point. This weighting function depends on the spacing 
of the nodes next to the point of interest and the 
interpolation function used in deriving the element. In the 
results presented the interpolation functions assume linear 


variation between nodes of all properties, including shear. 


‘For Newtonian fluids only. 
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Figure 5.1 and Figure 5.2 are for the case of a 
circular pipe to show the typical characteristics of the 
shear force F as K is varied. The pressure gradient is a 
constant. Figure 5.1 shows how the shear force on one node, 
at time t=0, varies with K. The value of F is proportional 
to the pressure gradient and the area of the pipe. The force 
per node is inversely proportional to the number of nodes 
used to represent the circular pipe. For example, if the 
area of the pipe is unity and a simple arc of 45° is used to 
model the circle(i.e., there are effectively 8 nodes on the 
circumference), the force on each node, at K=0, is 1/8 of 
the product of the pressure gradient and the pipe area. 

As K increases the force at t=0 drops quickly at K=2 
but at K=4 it seems to be decreasing asymptotically to zero. 
Figure 5.1 also shows the same example with wt=90°. As K 
approaches zero obviously the shear force approach zero. Its 
behavior is very interesting after K=2. It reaches a peak 
near the value of the force at t=0 at about K=2.5 and begins 
decreasing identically to the force at t=0. 

Figure 5.2 shows how the magnitude of the shear force 
|Fs| defined as 


LZ 


irepeee’ (rsr2*e<Poie) 5.2 


varies with: K. 


The mesh of Figure 5.3 is included to provide a better 


understanding of the following results. The symmetry of the 


He ie 


§ to peer Pee B | 
gas’ te cpixe Lute Hatin 2 
s wt ¢nsioatp re ee at ae | 

shoe 265 ae $2963 viigat it Ww ‘ ‘ic he ba enti : 
fenolitteqnig? et ¥ Ll ea j ape Soin ee 
doxoh att .agig add 34 eae wtts Bis dnotbert' 6 63 aoty 
aghon 30 vedere HEROOT Couey cet ibeaseen! at bel 
sc) i) , ake 204) aig i dare fa re aoa 
Ot bike "Ot Wo 938, elgmke Ban a ef stata ad: 
ett 10 eehad 6 (Lowitcgn fe onelvends a “ts tat 
- 30, el. ef. Cel ve abon soe at eset a teon ial ons 
ate ag ht a3 bere neha 239g ts qe ‘toate 1o% o 
Sah du yrkoiue aqoxb Cheesy, eae eC ; 


A eh oo Pope vty sonny + 


me 


rf Lis ae 
Ses ted @ eoive4g of banat gael wk be sxup it TY iene 2 
ait 6 yioomeys et? oo OGmD GReWel fot 243 to om 


| i‘. ine 


86 


triangle and the square are used to the fullest. 

Figure 5.4 to Figure 5.15 show the shear stress the 
fluid exerts along the pipes at various frequencies. Figure 
5.3 shows the mesh used to obtain the final results of the 
triangular and square pipes. It is the same mesh as that of 
Figure 4.21. The only difference in modelling the two types 
of configurations, as far as programing is concerned, is 
which sides of the pipe are constrained and the difference - 
of scale to make the pipe have the desired hydraulic radius. 
For example, to model the triangle of hydraulic radius 1 the . 
dimensions of the mesh in Figure 5.3 is used with the nodes 
labelled A to U constrained. To model the square of 
hydraulic radius 1 the scale of Figure 5.3 is reduced by 
2.414 and only the nodes Q to U are constrained. 

In Figures 5.4 to 5.15 each data point on the graphs 
represents an average Shear stress inethe area of a node. 
The data points are labelled to show data points 
corresponding to each node in Figure 5.3. 

For Figures 5.4 to 5.15 the vertical axis is shear 
stress for the pipes with hydraulic radius of unity and 
pressure gradient of unity. For Figures 5.4 to 5.9 the 
horizontal axis is the distance along the perimeter of the 
pipe beginning at the corner of the square pipe and ending 
at habf athe wwidtha:Seeapoints OrctocUtel ehuquressa26 For 
Figureso5.106to 5.15 the horizontal axisers; the: distance 
along the perimeter of the isoceles triangle beginning at 


the right angle corner and ending at the midpoint of the 
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opposite side of the triangle. The distance from A to U is 
the sum of A to Q and Q to U. Again, by symmetry, 
information of the variable is completely shown by these 
cases. 

The case of K=0 1S worth some extra comment. When the 


frequency is zero the governing momentum equation reduces to 


which is of the same form as Prandtl's torsion function for 
prismatic bars from the field of elasticity. (see 
Sokolnikoff [17] or Timoshenko [18] ) It can be shown that 
the shear stresses of the steady state fluid flow are 
analogous to the stresses of a bar in torsion. The largest 
stresses occur at the perimeter of the configuration in both 
cases. 

It may appear as an obvious flaw that in these figures, 
at K=0, the shear stress is not zero in the corners. Since 
the matrix equation produces the shear force in the area of 
the node and dividing by the area represented by the node 
only gives average shear stress it is not possible to define 
the stress in the corner. A finer grid along the perimeter 
of the pipe will produce more accurate shear stress results. 

From these results, however, some general observations 
may be made. At low frequencies the shear stress is nearly 


parabolic along the flat edges of both the triangular and 
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Square pipes. As the frequency increases the profiles 
flatten and the magnitude decreases. In all cases the stress 
approaches zero in the corners of the pipes. 

The primary use of the shear force matrix equation is 
to find the average shear stress on the pipe in order to 
calculate the Biot curves. The average shear stress is found 
by canta the shear forces supplied by the matrix equation 
and dividing by the length of perimeter they act on. This 
procedure should be more accurate than eaten tens a number 
of shear stresses around the pipe and averaging them because 
of errors introduced in finding derivatives of the velocity 


numerically. 


‘Ba 


‘viel Fresg ddl eveyome 4 a iv hep 
tate ets eoana Dia at re seed 
weqig ace to. eames) mi exe as 

2i eolssupe «i¢cem era 4 apni: i ae SR. seairy a 2 Ai 


_ 


o3 2ab20o af Soi @ SR Re sae any ne penser. oe aisd 
= ie 


brval ef 2am? fe, came aeene"s ant eeovaus ‘Soke ~ad¢ sg 
eaitcuns xlacanrets 2 bri lege auarar, ‘taeds ott vie | 
eayiT «no fon Waly resem Ney *” rote wis onset : 
redi@n, 6.206 ahs shy ned? etde foo he OMA at Bluede ou 5930! 


Hmecsa wand go er Uern arts aie eds iauoaae etaeavie. 4) | 
gtionley ody leo eal gy 5 eal ‘gee ni wiaciaibs wi0% oat 


i‘ a4 i 
1 ‘ 
{ wi 

"4 iu F 
Ae | pee 

7% 
{ quinn “ee i Lge 

| fv, ; [ f 

eo.) / 


j 7 «af 
boost 
& e 5 ; _ ? 
Ne t » 7 
| ap mu : F 
Le , 2 DIA ee 
Jie } a 
ae Ni 
ue : 
) re AC he 
, i 7 an hm ; cA) 
5 ' 
"i . F, 
- ‘i 7 
al eae 
. ’ ora ; 
a) 2 ' 


-0.0025 0.0 


-0 »0100 


o—o shear at t=0 


o—® shear at wt=90° 


SHEAR @dP/0x ) 


-0.0175 


v4 


VA 


-0.0250 


0 1 2 3 & cs) 6 a 8 9 10 
K 


Figure 5.1 Shear vs. K for circular pipe at two times in 
the cycle 
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vs. K for circular pipe. 


Figure 5.2 Magnitude of shear force. 


Z ~ oy ity) F 


ie 4 u 
ou : 


7 e454 fl ALY = . 
; 1 Ae ot ee ) 7 * a y ’ 
an aes es Sk 
ileal, oo arbietan tied cnn 
4 er Te? en r* wo 
eee”, | ' : 


91 


PONG DE Se eye. El. 


ici, 


Figure 5.3 79 node grid to model triangle and square. 
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Figure 5.4 Shear stress for square pipe.Lettered values 
correspond to nodes of Figure 5.3. 
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Figure 5.5 Shear stress for square pipe. Lettered values 


correspond to nodes of Figure 5.3. 
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Figure 5.6 Shear stress for square pipe. Lettered values 
correspond to nodes of Figure 5.3. 
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Figure 5.7 Shear stress for square pipe. Lettered values 
correspond to nodes of Figure 5.3. 
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Figure 5.8 Shear stress for square pipe. Lettered values 
correspond to nodes of Figure 5.3. 
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Figure 5.9 Shear stress for square pipe. Lettered values 
correspond to nodes of Figure 5.3. 
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Figure 5.10 Shear stress for triangular pipe. Values 
correspond to lettered nodes of Figure 5.3. 
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Figure 5.11 Shear stress for triangular pipe. Values 
correspond to lettered nodes of Figure 5.3. 
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Figure 5.12 Shear stress for triangular pipe. Values 
correspond to lettered nodes of Figure 5.3. 
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Figure 5.13 Shear stress for triangular pipe. Values 
correspond to lettered nodes of Figure 5.3. 
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Figure 5.14 Shear stress for triangular pipe. Values 
correspond to lettered nodes of Figure 5.3. 
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Figure 5.15 Shear stress for triangular pipe. Values 
correspond to lettered nodes of Figure 5.3. 
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6. Pipe Acoustic Properties 


6.1 Resistivity and Effective Density 


This analysis began with a simplification of the 


Navier-Stokes equations 


du #) 
X eo. Fo eH Vee. 


from which are found the velocity distributions in the pipe. 


Now uSing only the mean velocities 


= | XE 
Ue ads dare 5 (Uh da) 6.2 
and writing the above force balance as 


a= (5 w +R) U 


Q 


where R is the resistivity per unit length (Rayls/m in MKS), 


Peis the effective density. Rearranging the above equation 


results in 
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It is known from flow dynamics that for steady state 


flow (K=0) ina elie slit of width b 
Re=o20 b+ | 6.7 
and in a circular pipe of radius a 
R = 8 /a? “Wig G58 


where u is the dynamic viscosity. The results of the model 
match these values very well. 

This information is very important to the application 
of the results. To determine the resistivity for any pipe at 
any frequency the following steps should be taken. Determine 
which of the example shapes best approximates the shape of 
the new pipe. Calculate the hydraulic radius of the new 
pipe. Calculate the K from the hydraulic radius, the fluid 
properties and the circular frequency of interest. At this K 
find from Figure 6.1 the R value. Multiply it by the dynamic 
viscosity and divide by the square of the hydraulic radius. 

With a slight variation of the above procedure it might 


be within reason to find approximately the resistivity for a 
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porous material. The following has not been verified, 
however. The first Heo we to experimentally find the 
resistivity of a sample of the material. This is done by 
finding the pressure drop across a sample of material (in a 
pipe where the sample seals a portion of the pipe) given a 
known steady flow of air through the sample. This 
information is used in equation 6.5 . Since the frequency is 


zero 
R=- 1/U(dP/dx) a7 


Using Lord Rayleigh's [15] assumption that a porous material 
can be approximated by a matrix of parallel channels then 
from equation 6.7 a typical pore dimension can be 
determined. With this information the K can be determined at 
whatever frequency is of interest and R can be scaled as 
described above. 

As mentioned the extension to porous materials needs 
testing. A simple test with a few samples could be done to 
see if there is good correlation to pipe flow. One 
difference could already be predicted and that is the 
position of the peaks of absorption if an impedence tube is 
used. The path length in a porous material is greater than 
simply its width so the peaks would be shifted to the lower 
frequencies. Nonetheless the procedure would be very useful. 

For very high frequencies Zwikker and Kosten [19] have 


shown that 
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The finite element results approach these values for the 
Circular pipe. 

Crandall [3] has shown that as the frequency approaches 
zero, the effective density in a circular pipe approaches 
four thirds the static density. Again this is the behavior 
of the numerical results. As the frequency increases the 
effective density approaches static density for all 
geometries of pipe. The value of the density is dependent 
only on K. Unlike the resistivity no scaling is needed. 

Figure 6.1 shows the resistivity, R, for various 
geometries of pipe. The viscosity used is unity. In each 
case the hydraulic radius is also unity. Figure 6.5 shows 
the relative size of these configurations, each with 
hydraulic radius of wnacty. 

Figure 6.2 shows the ratio of effective density to 
static density under the same conditions as noted above. 

To note the sensitivity to the grid type Figure 6.3 and 
Figure 6.4 are included showing the results of the same 
variables. The rectangle's grid in Figure 6.3 and Figure 6.4 
is uniform whereas in Figure 6.1 and Figure 6.2 the 
exponential transformation, similar to equation 4.4, is 
applied to concentrate the nodes near the walls and corners 
of the pipe. In Figure 6.1 and Figure 6.2 the grids for the 


trangular and square pipes have approximately twice as many 


rt 


ho GBs 


uli? a 


frat svioding ae ots: wert ea. 
an bd choo) sith ine an, aebis caienetl: > as - 


¥ ua 
Rn 4 28 270 


nate p44 


og ».9 


ie. 


560705 DOs 


ae EO Bsr. wi ea 


iy shades antanartnted Ee 


mnie 32 pdtv dal, Shes Leingiias vies PH 


eydsupetT eas a oT AY: ‘SPs: ero 


‘nad oh Tee of eaiiallar ti oes 
? tay BP: Sagu cite iro ant vaciq ie @ 
vive lS «824 orbs at ‘hula: te 
‘Cha¢-q nal sa LRH. thats Ae, waka 


e Seen ay ok 
‘ my A rm on a | 


i 
J 


é c ae 


x 
: : ~“ 
Dae my as oh 


a ee setae beipsettw a 
S144 °rivod ‘SAKES. 9 vere y od 88 ; 


diene! Reese fetsausia es artim 


giovrsh: eit I syne oul‘ ee 34 lay 
arilease wn 9 vedanta oooh ap ab aity wht 


engand rt au tbee 4 th | r 


Te) : 


{ue 


agg ih — haste Tienes mt + 08 


ar pay ri nage 
soupe ot Sahih pie: fearon) pti: 
z pees 8 sramtanitie a ® 
Pet s3Upey A uepba #2 30) 


fie met 


Sree mi a) +9 i: Bigic ‘erenpe baw hpatmie as a 


fa? Ra “W tap ag 


108 


elements. In both eas the type of spacing of the nodes is 
the same i.e. the elements are relatively finer to the same 
degree near the edges. 

This indicates, first of all, that the density term is 
not greatly affected by varying the grid. Second, and more 
importantly, the resistivity is Significantly raised at high 
frequencies by using a better grid. Thus it should be noted 
that even though from Figure 6.1 it seems there is 
Satisfactory convergence the program does yield a lower 
bound to the actual resistivity. 

It 1S interesting to note the ordering of the curves. 
The order of the resistivity curves is reversed compared to 
the effective density curves. In addition the order is 


unrelated to the order of the Biot curves. 


6.2 Impedence via Exact Equation 

In order to completely specify the acoustic 
characteristics of a material two independent functions are 
required. In the previous section these were the resistivity 
and effective density. A more common method to describe the 
characteristics is with the complex impedance Z. 

The impedance is based on the particle velocity 
generated by a given sound pressure at the surface of the 


absorbing material and is defined by 
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In this case the absorbing material is the matrix of 
parallel pipes and the "Surface" is at the open end of the 
Pipes. The other end is assumed to have a rigid terminus to 
Simulate the specimen being mounted against a hard wall. 
Appendix A derives the relation between the impedence 


and the absorption coefficient a as being 
a = 4Real(Z) / (|Z|? + 2Real(Z) + 1) 6.13 


Lines of constant absorption are anbun in Figure 6.6 . When 

Zr iS unity and Zi is zero complete absorption occurs. 
Appendix A also contains the derivation of the 

impedence in terms of the resistivity and effective density. 


The result is 


Z = Zo coth(jBpL) 6.14 
pce 


0 
e : L/2 
ee Satay fee aN RS) 
) 8B pees j R/o.w) 6 
ce) 
This is the exact answer assuming the pipe is of constant 
cross-section along its length. Figure 6.7 shows the 
impedence curve using the above exact equation for an 


isoceles right triangle with hydraulic radius of 0.00023 m. 
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6.3 Impedence via Acoustic Finite Element 

A second finite element program was developed which 
approximates the pressure distribution along a pipe given 
the resistivity, effective density, the frequency and the 
input velocity. 

This is done in such a way as to give us the ability to 
model the pressure distribution in a pipe which varys in 
area along its axis. This is of some advantage over using 
the exact equation of the previous section because a pipe of 
varying | area may more closely model a porous material. This 
is suggested for further work. 

The derivation of this program is based on Galerkin's 
method also. Craggs [2] derives a similar element from 
variational principles. 

The initial equations are the force balance and 


continuity, respectively 


- grad P = (j pw + R) U 6.16 
- 7 P = div U 6.17 
ee 


In matrix form the final equation is 


@) ; 
2 e R 
[s] - kK a (Py - k 5c [P] LP} {9} 
, ‘ Bey 
0 Por 
BR i ke Je Pp a 
ao eel [s] - k mt {P} {Q. 
{e) 
6.18 


where {Qr} is found from R, the interpolation functions and 
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the input velocity at the open end of the Pipe. The array 
{Qi} is similar but in the place of R is the product of the 
effective density and circular frequency. [P] is referred to 
as the acoustic mass matrix and [S] is the acoustic 
stiffness matrix. 

Appendix D contains the derivation of the above 
matrices and matrix equation. At first [P] and [S] are 
derived using only the above equations. However, the varying 
area can be taken into the calculation of these matrices by 


introducing the horn wave equation, which is 


iy etig vt 3 | 9P 
a eye) 3z {3x M00] 6.19 
se C2 A(x) 0X [dX 
This equation leads to 
it 
mh nay Cue 
[s] = f og Screeo. ax 6220 
@) 
[el =) eeneeentD 0) dx 6.21 


where L is the element length. 
This program was tested with interpolation functions 
that assume linear variations of the area and pressure 


between the ends of the elements, for example 


Pixie P61. - x/l) +e Peek bh 6.22 


P, is the pressure at the node at x=0 and P, is the pressure 


at the node at x=L. 
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These interpolation functions resulted in 


[s] = =e i | 6.23 


It is suggested that, since the advantage of this 
program is in modelling porous materials a higher order 
Pe tion should be used thus eliminating at least the 
discontinuities of the derivatives of the area. For now it 
is sufficient to use the linearly varying acoustic elements 
Since the exact equation (used for verification) is only 


valid for constant area pipe system. 


6.4 Results 

The results of the acoustic finite element model are 
compared to the results from the exact equation in Figure 
6.7 . In both cases the calculations are based on the case 
Of a, triangular pipe, of hydraulic radius of 9.23 mm. , 

The prediction of the model follows the trend of the 
exact equation less than expected. Craggs [2] derived a 
Similar element but obtained better correlation. 

The acoustic finite element (AFE) results are all 
closer to the point (Zr=1,2i=0). (at which complete 


absorption occurs, see Fig. 6.6) than the exact equation. 
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Comparing Figure 6.7 Figure 6.8 shows that doubling the 
number of elements does not significantly improve the 
results. The solution appears to have converged. Better 
results were expected because Craggs [2] with a Similar 
element (in a Slightly different use) obtained good results. 
Thus the predicted absorption is greater than that of the 


exact equation. 
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Resistivity (R) 


SQUARE (79 nodes) 
CIRCLE 
(79 nodes) 


i 
RECS (fine mesh) 


Figure 6.1 Resistivity for various geometries using 
best meshes. Hydraulic radius is unity 


for each case. 
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Figure 6.2 Effective density for various geometries using 


best meshes. 
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SQUARE (49 node) 
CIRCLE 


TRI. = 
SLIT (68 node) 


RECS (uniform mesh) 


Resistivity (R) 


Figure 6.3 Resistivity for various geometries using second best 
meshes. Hydraulic radius is unity for each case. 
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Figure 6.4 Effective density for various meshes using 


second best meshes. 
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Figure 6.5 Pipe cross-sections scaledto same 
hydraulic radius. 
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Figure ©-6 Lines of Constant a 


149 


120 


Siang PENT TEDELER: 


©——® EXACT IMPEDENCE 


4000 Hz 


Figure 6.7 Impedence for triangular pipe 
(Length = 0.076 n, r= 0.0023m) 
Calculated via exact equation and acoustic 
finite elements (10 elements). 
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Figure 6.8 Impedence for triangular pipe 
(Length = 0.076 m, r, = 0.0023m) 
Calculated via exact equation and acoustic 
finite elements (20 elements). 
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7. Absorption, Theoretical and Experimental 
Most of the results obtained by the computer model are not 
easily verified by experiment, eg. the velocity distribution 
in a non-circular pipe. There is fortunately a simple, 
easily repeatable experiment with which to find the 
impedance and absorption for a porous absorbing material. 


That is the impedance tube experiment. 


7.1 Absorption from the Numerical Model 

Our numerical model so far has led to the specific 
impedance, Z , of a pipe given its length and dimensions of 
its cross-section. From the specific impedance can be found 


the absorption (see also Appendix A) via 
auec4eRea ize yn | Zii79t 2Real(Z) + 1jpai 


These results can be compared to the experimental 


measurement of absorption by a mesh of triangular pipes. 


7.2 The Experiment 

The apparatus involved in the impedance test is seen in 
Fig. 7.1(a) and drawn schematically in Figure 7.1(b) . 

A plane sound wave, generated by the loud speaker, is 
reflected by the specimen at the opposite end of the tube. A 
partial standing wave results in the tube. A microphone 


attached to a moveable probe finds the maximum 
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and minimum pressures of the Standing wave. From this the 
absorption can be found, as derived by Kuttruff[12] (and 


summarized in Appendix A). The absorption a is 


a = 4(PmaxPmin) / (Pmax+Pmin) ? 7.2 


The distance between the specimen's surface and the 
first pressure minimum is required to find the phase angle 
in the reflection factor. Denote this distance by Xmin. The 


phase angle is 


where A is the wavelength. 


The reflection factor Rf is defined as 
Re Ma IR exp (jw8) 7.4 
The absolute value of the reflection factor can be 
found from the maximum and minimum pressure also 


| R£|= (Pmax - Pmin) / (Pmax + Pmin) 7-5 


To find the specific impedence 
PesesiteRe ) o/) (1 -oREt) i26 


The dimensions of the impedence tube fix limits on the 


range of frequencies which can be reliably tested. The 
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length of the tube must be long enough to permit the 
formation of at least one maximum and one minimum of the 
pressure distribution at the lowest frequency of interest. 
For the following results a tube of approximately 1 meter 
was used. This set the lower limit around 200 Hz. Actually 
the results between 0 and 400 Hz are interpolated (at f=0 
attenuation is zero). The diameter of the tube set an upper 
limit on the frequency. For circular tubes the diameter must 
be less than 0.586 times the wavelength (see Kuttruff [12]). 
In obtaining the following results a tube with a diameter of 
0.05 m was used. This set the upper limit at 4000 Hz. 
Kuttruff [12] mentions that since loudspeakers are not 
generally free from distortion and therefore higher | 
harmonics of the meaSuring frequency are produced, causing 
errors, it is recommended to use octave filters to remove 


the frequencies not desired. This advice was followed. 


7.3 The Specimen 

The sample used to obtain the experimental results was 
provided by General Motors Corp. 

The specimen is made by laying a sheet of finely 
corrugated foil on a layer of plane foil and rolling up the 
two layers. The result is essentially a matrix of parallel, 
nearly triangular, tubes. See Figures 7.2(a), 7.2(b) and 
7.3 . These pipes are described as "nearly triangular” 


because it is not possible to bend the foil in precisely 90° 


> yt 


Marit oe 


ie eee se sila, aint 


my 2) €2534R- 8 Sauna Se eae — aaah ile 


Dali! — 


ae 

Gs - " : ; 
am! ' Yee, aw iy oly 
| 7 7 =. Ha 
an ¥ ee i Q an | <? 


ee 


ua hs ban, hee eat na 8 ete sqedi, wg 
8 eg niga? See we Ss ee “ha onbd Suites 

an. eorqe am gus r sghunigsoniaek BOR <a 
APA us ean oa4 soa gaeNe: 2 

eosnd. ou it Bos tae amgwndcses cult 


— '- 7 
a0 me = 


—. 


e: 507 Oa EB aee caer, ee ag tany 4 
aseud . caller i ae Crhengens oak Ae ryt 


04 AS9nGs tv 5, ORS bore ane. * ay azaty 


ati egy e4o.296. eca?e . bet “av ma 


yon afi asi 1se ie! 4 3 tue oy 
tiv 29560949 U9 RELIC Jee re $3'r3 nts a 
g° wtsot pol ene, att So: mabe Taub 
; 4 


ene 
yaverée eh, ra ive ight os oe $y ee 


— § i . 
naw” eos cba ghad Whesi 0k on ee ag 


i | a 


oe Mu és a | =) 
et Vs if ee “omen att | Le 
» fagnaat bon ea bogie nna ven = 


shi So Sawae inet or ehaia @) at hace AR. yan 
aarti l+s, Ofen tint ite te te9el #0 ceae bebact mS 


© 
va uy 


125 


corners. There will always be a non-zero radius of curvature 
at the corners. It is assumed that a typical pipe in the 
specimen can be modelled by an isoceles right triangle with 
a hypotenuse of 1.1 mm. This means its hydraulic radius is 
0.23 mm. The specimen of Figure 7.2(a) is 76 mm in length 
and 51 mm in diameter. 

For oe Sample in air the 4000 Hz limit of the test 


equipment translates to a K of 9.2 in the triangular pipes. 


7.4 Results 

General Motors Corp. also supplied their own test 
results for the specimen. Figure 7.4 compares present test 
results with those supplied. Since the curves are very close 
to each other it 1S reasonable to conclude that the 
repeatability of the impedance tube method is good. In the 
remaining figures the experimental data is from our test. 

Figure 7.5 and Figure 7.6 compare the experimental 
results to the absorption obtained using the impedance of 
the exact equation. The only difference between the results 
of the model is the speed of sound used. The values used are 
C=290 m/s and C=343 m/s respectively. 

The experimental and predicted results compare very 
well. The slight discrepencies between the results can be 
broken into two categories; the height difference of the 
curves and the frequency at which the peak absorption occur. 


The height difference depends on the length and area of the 
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Pipe, the limitations of the assumtions made on the fluid 
behavior and the limitations of the numerical calculations. 
The frequency at which the peaks occur depends on the length 
of the pipes and the speed of sound in the pipes. 

To model the corrugated foil specimen it was assumed 
that the pipes were similar and could be represented by an 
isoceles triangle with a hydraulic radius of 0.23 mm. A 
small error in meaSuring the hydraulic radius of the pipe 
would cause an error in the absorption of the same order of 
magnitude. 

In deriving the model of the fluid motion an assumption 
1s made in neglecting the temperature changes so as to allow 
the use of a constant viscosity. The assumption is that the 
wall of the pipe maintains a constant temperature and the 
heat conductivity of the fluid to the wall is perfect, this 
means isothermal compression and expansion. It can be 
reasoned that this is closest to the truth at low values of 
K. At higher values of K the process becomes adiabatic. 
Increasing K can be thought of as increasing the radius of 
the pipe, this means the fluid in the center of the pipe can 
less easily transfer its heat to the wall. With less heat 
transfer the process becomes adiabatic. 

If the compression is isothermal the speed of sound in 
the pipe is C=290 m/s. The first absorption peak of the 
numerical model for C=290 m/s is at the same frequency as 
the experimental data. It can be shown that the frequency at 


which the peak occurs is proportional to the ratio of the 
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speed of sound in the pipe and the length of the pipe. 
Therefore since only the first peak matches with the 
experimental result this is confirmation that C=290 m/s at 
low frequencies and the compression is isothermal as 
assumed. The second peak of the numerical result in Figure 
7.5 is left of the experimental peak. This move to the right 
of the experimental results means the speed of sound is 
increasing. Figure 7.6 shows the effect of using C=343 m/s 
fon the model. 

It 1S possible to improve the matching of the peaks of 
the absorption curve by using empirical relations to 
determine the speed of sound at various frequencies. For 
this sample this would become more important at frequencies 
above about 3000 Hz. 

Bamb oie, in approof which, is aivariation of that 
given by Lord Rayleigh [15], derives the decay of a wave in 
a circular pipe considering viscous effects only. He 
mentions that a more complete investigation was provided by 
Kirchhoff [11] in which thermal processes are considered, as 
well as viscosity. Lamb [13] compares the different results 
of decay and states that the decay is increased by 
considering thermal processes but remains of the same order 
of magnitude. 

No consideration is made in this model to account for 
dissipation due to thermal effects or the viscosity changes 
G@tcher across therpipe or with respect, to thercycle. ‘These 


considerations are neglected by the assumption of constant 
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temperature at all time. This assumption is likely the 
largest single factor responsible for the deviation of the 
absorption of the model from experimental results. The 
assumption, as mentioned, is best at low K as evidenced by 
the experimental absorption increasing faster with frequency 
than the absorption of the model. 

Another factor to consider is the limitation of the 
numerical technique. When comparing the resistivity curves 
of a particular pipe modelled with different grids, in 
chapter 6, it was noted that the model produces a lower 
bound on the resistivity. In addition the accuracy decreased 
with frequency producing a much lower bound in the case of a 
crude mesh. Thus the calculated absorption will be a lower 
bound and will be further from the experimental data at 
higher frequencies. 

At higher frequencies there are two factors causing 
drift in the same direction, the lack of thermal effects and 
the limitations of the mesh. A number of meshes were 
compared and the one used to generate the final data seemed 
to be very close to a converged answer. Therefore it is felt 
that most of the error, albeit small, is due to the lack of 
consideration of thermal effects. 

A point of discussion on the absorption obtained from 
the model is why the height of the absorption curve is 
greater for C=290 m/s than for C=343 m/s. From the theory, 
the only influence the speed of sound has is on the location 


of the peaks, rather than the value of the absorption. 
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In conclusion,from these comparisons with the results 
of the experiment it seems the integrity of the analysis is 


very good. 
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(a) Testing Equipment. 
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(b) Schematic of Testing Equipment 
Figure 7.1 Testing Equipment 
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8. Discussion and Summary 


8.1 Method of Solution 

In defining the problem, necessary assumptions are made 
to reduce the problem to a manageable form. The Finite 
Element Method was chosen to solve the defined problem. This 
method is relatively new to the solution of fluid motion 
problems. Most of its use is in the field of solid 
mechanics. The terms ‘stiffness matrix' and 'mass matrix' 
Originating in solid mechanics have been adopted here 
because of the analogous use and behavior of the fluid 
motion to structural motion. 

The method of solution initially had the requirement 
that it must be able to handle any shape of domain. It was 
later found that for efficiency it also must possess the 
flexibility of concentrating on the solution near the edges 
of the domain. The Finite Element Method is thus 
particularly applicable to this problem. Fairly good results 
can be easily obtained on existing computing systems. 

In defining the problem an assumption is made about the 
form of the pressure wave being harmonic. This assumption is 
of convenience, not necessity. Even with this assumption the 
solution for an arbitrary input pressure function could be 
found by adding the solutions of its component frequencies 
and amplitudes similar to Fourier series of wave forms 


because the governing equation is linear. If for some reason 
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this procedure is not desired it is valid to begin with the 
differential equation (Equation 2.2) and define 
interpolation functions with a dimension in time i.e., 
Ni=Ni(y,z,t) and input nodal pressure gradients as a 
function of time. This might be the best way to find the 
velocity distribution due to a pressure gradient which is, 
for er anpies triangular with time. 

Galerkin's method is used to derive the matrix 
equations. This method of approximation is powerful, 
relatively simple and can be applied to any linear or 
non-linear governing equation. The equations resulting from 
the Galerkin's method are the same as those from the 
variational principle. 

Galerkin's method was chosen first to benefit from its 
advantages but also in order to become familiar with a 
technique which is readily applicable to many problems. 

This analysis produces an approximate solution to the 
velocity profile in any shaped pipe at any reasonable 
frequency and at any time in the cycle of the periodic 
pressure gradient. These results can be directly and easily 
used to find the shear forces on the pipe wall, the stress 
State of a bar in torsion, as well as the acoustic 


properties. These properties include the resistivity and 


'Galerkin's method would produce different equations, though 
still valid, from Lagrange's equations if the original 
equation had been multiplied throughout by an arbitrary 
factor. However, in mechanical problems there is an optimum 
way of applying Galerkin's method which renders it 
equivalent to the use of Lagrange's equation(See Ref 6). 
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effective density, and, given the pipe length, the impedance 


and absortion. 


8.2 Conclusions 


1) 


2) 


3) 


oy, 


6 ) 


The major points to be noted are: 

The principle of similarity can be applied to the 
PesiStivityeresults) so that the acoustic: properties (of 
most configurations and any size of pipe can be 
determined @iron them iniormationy presented. The ertective 
density results do not require scaling. 

The correlation of the experimental data to the finite 
element data is very good. 

For the cases of resistivity and absorption the program 
gives a lower bound to the solution. This is most 
relevant when uSing results from a relatively crude 
mesh. 

At K=0, the solution is applicable to finding the stress 
state in prismatic bars subjected to torsion. 

Accuracy is most dependent on concentration of nodes near 
the pipe wall. 

The Acoustic Finite Element can be used to model 
impedance of a pipe with varying area along the axis. 

The analysis and method of solution have been verified by 
experiment. Problems with similar assumptions, governing 
equation and boundary conditions can be solved with 


confidence via this method. 
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8) Some improvement could be made in the model by accounting 
for the influence of thermal processes in the decay of 


the wave. 
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Appendix A, Acoustic Properties Equations 


This section will show the steps for deriving the 
absorption of a wall in terms of the specific impedance and 
in terms of the maximum and minimum pressures(to be used for 
the impedance tube experiment). 

Assume that the wall is at x=0, the wave is travelling 
parallel to the X axis and the wave direction is normal to 


the wall. If the incident wave is coming from the negative X 


direction its sound pressure is 


Pi(x,t)=Po.expl[ j(wt-kx) ] A.l 


For a plane wave moving in the positive X direction 


oP i dU 3 A.2 
Therefore, the particle velocity of the wave is given by 


= ) . 
a wear j -k A.3 
Vert) 6c exp{j (wt-kX) } 


For the reflected waves 


Fete) wees exp{j(wt + kx) } A.4 


where Rf is the reflection factor, 


<—_ Po 
-R am exp{j(wt + kx) } A 
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Re = IR, exp (38) A.6 


Both |Rf| and ® can be found from the impedance tube as 
described near the end of this section. The reflection 
factor 1S a property of the wall and describes both the 
attenuation and phase change of the reflected wave. 

The total sound pressure and particle velocity in the 
plane of the wall are obtained by simply adding the 
expressions for the incident and the reflected waves and 


setting x=0 


P =p +R 5 : 
(6.4) >t f) exp (jwt A.7 
a 
eve = ae ~ R,) exp (jwt) Awe 


P(0,t) divided by u(0,t) gives the impedance 


1+R 
oaeuzes (o,t) _ pC oo Lo 
. (0, t) ° 
Conversely 
oe A.10 
Eris iz ca 


where Z is the specific acoustic impedance. 
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Since the intensity of a wave is proportional to the 
square of the pressure amplitude and the absorption 
coefficient a is the loss of incident intensity on 


reflection it follows that 


2 
a=1 - |R,| Ratt 


Substituting eqnAl0 into eqnAii1 gives 


Re (Z 
ane e (Z) INA ILA 


\z|7+ 2 Re(z)41 


The relation between the absorption coefficent and the 
specific impedance can be seen in Fig.6-6. 

The distribution of the sound pressure in the standing 
wave due to the wall is found by adding eqn 1 and eqn 4 and 


taking the absolute value 


‘ 2 3 1/20 
Phir rates IR, + 2|R,| eos Ce A.13 


The pressure amplitude varies periodically between 
| Pmax=Po(1+|RE£ |) A.14 
and 
Pmin=Po(1-|Rf |) A.15 
The distance between maxima is half the wavelength. 


From equationsA11,A14 andaAl5, 
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2 max min 
IRe| po Pp K.16 
max min 
“ ae fF Pears 
Qq = 5 IN SALT 
(P + Pp ) 
max min 


In the impedance tube it is easy to measure the minimum 
and maximum pressures as well as the distance between the 
surface of the test specimen and the first pressure minimum. 
This is all the informatiom necessary to determine the wall 
impedance. 

From eqn 13 the first minimum occurs when 
cos(2kx+r)=-1. Denoting the minimum distance which satisfies 
this condition by Xmin the phase angle 8 (of eqnA6) is found 


from 


where lamdaiis the wavelength=c/(frequency) 
Thus the absorption, reflection factor and the impedence of 
the sample can be obtained from the impedence tube. The 
following, from Craggs [2], shows how the impedence can be 
derived in terms of the resistivity and effective density. 
The exact solution for the input impedance, 2, of a 
column of absorption material of length 1 with a rigid 
Beeneeicn is given by Zwikker and Kosten(9) as 


Z=Zo.coth(jlk), where Zo is the characteristic impedance and 
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k is the complex propagation constant. The quantities Zo and 
k may be obtained considering the one dimensional form of 
the equilibrium and continuity equations, as shown by 
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plane harmonic progressive wave 


“ a 


jKP - (jup, +R) U=0 rival 
ay “—s \ap) 4K Ut 70 A.22 
pc ‘ 


For a non-trivial solution the determinent must be zero. 


A.25 


jK -(4 + 
J (Jwo . R) 
nO) A.23 
ee ; jK 
€ 
Po 
This gives 
p F 
W e jR 172 
is ye araea lin See) A.24 
(es e 
ae Pion 
d 
an x & 
ae 1/2 


oa) | asin eon tn 


_ ae 


a i 


S 


me. - = re 
. ’ ’ ie si a ae ——/ 


oles ht bab siiataewale 
bide \ditbineutety ssiclealalt 


\ ~ 
i i et! i ae ae 
a ‘ 
a | NS » gli 
a a ; 
wi er 
} 
ae 
; a nl = 
a 
Oy 7 
oy -~ cf at 2 7 —— {* Tt) i) 8 
7 T 
é ; 
Pi | a i 
a a) 
: 7 / 
wi ad 


147 


Appendix B, Fluid Finite Elements 


This appendix is to provide the discussion of the steps 
necessary to reproduce the finite element solutions or 
provide a basis on which to build extensions of the 
analysis. 


The momentum balance governing equation is 


From the wave equation we find the satisfactory form of a 
harmonic wave 
A 
P=P.,exp(j(wt-kx) ) Baz 
Substituting a velocity of this same form into eqn B1 gives 
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av) 
Let us assume U to be in the form 


Ni are interpolation functions, Ci are the unknown field 
variables to be found. Ci are velocities at the r points 
in the domain if linear interpolation functions are used. 

If the order of the interpolation functions is 
increased, Ci must include derivatives of Ehemiield 
variable. In the following solution let Ci=Ui. 

To make the residual small in some sense we may set the 


total error’ to zero, Lue. 


figs dA = 0 Bao 


But to solve for the r unknowns we use r independent 


weighting functions in the above integral, i.e. 
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Galerkin's criterion specifies that the weighting functions 
should be the same as the interpolation functions. 
Therefore, applying this to our governing equation (B2) 
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Integrating term by term we define {P} and [M] by 


(2 Pda Py B.8 
pw fNPLNI aa (U} = pw[] {ov} ee 


The method of weighted residuals requires that the 
order of the interpolation functions must be the same as the 
highest order in the equation to be solved. Thus to reduce 
the order of our equation we integrate the Laplacian on U by 
parts. This, in most cases, conveniently introduces the 
boundary conditions as well as reduces the required order of 
the interpolation so that linear variation between unknowns 


may be assumed. 
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re) ) 3 oN- 
uJ tl ot pao a + foe pee) aa fu} = ull tu} Bet 


Now in the condensed notation we may write 
{P} =.(u(D] + Jp ow [M]) fu} Be12 


and taking on the notation of chapter 2 
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The solution of U now involves solving the above 
complex matrix equation. The variables {P}, {F} and {Uu} are 
complex as well as the term in the round brackets. 

Some computing facilities cannot easily solve for 
variables comprised of real and imaginary values. In order 
to solve for variables with only real parts we take the 


example of solving for U in the above equation 
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Breaking the variables into 
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Multiplying term by term and separating the real terms into 


one equation, the imaginary terms into another gives, 


respectively 
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This is the equivalent of 
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and we have eliminated j from the equations. Treating the 


shear force portion of equation (B13) similarly gives 
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chosen interpolation functions assume linear variation of 
the dependent variable between nodes. This means that the 
only unknowns at each node is the velocity i.e. 
consideration of derivatives is not required. 

Generalized coordinates were tried but this system does 
not work for some orientations of the elements. This problem 
waS circumvented by adopting a natural coordinate 
Syeten_ see Huebner | 2)¥ ror derails. 

On comparing the results of the two methods it was 
noticed that for most elements the elemental matrices 
produced by one method are entirely different from those of 
the second method. Only for elements that are square in the 
YZ plane are the elemental matrices identical. In spite of 
the different elemental results the two methods produced 


identical velocity profiles. 
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Appendix C, Biot Equation 


The following is a summary of the work done by M.A. 
Biota (i). 

Biot [1] Hemel the case of two parallel plates as 
well as a pipe of circular crosssection. This summary 
contains only the case of a circular pipe. 

Consider only the two-dimensional motion and neglect 
all pressure gradients and velocity components normal to the 
boundaries. The x direction is parallel to the boundaries 
and the r direction is normal to it with the boundaries 
being r=a. Let the walls be rigid and stationary. The xX 
component of the velocity is u(r). 


The equation of motion is 


where p is the fluid density and Vv’ is the Laplacian 
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factore?** , RewritingC-3 without this factor, 
BP ate eeiiooe 0, du ee 
OX JP E 52 ay j 


This is Bessels differential equation for U. The solution 
for U, remembering U is finite at r=0 and the boundary 


condition U=0 at_r=a, is 
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deviation from Poiseulle flow friction as a function of the 


parameter K. For large values of K 
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Appendix D, Acoustic Finite Elements 


It is implicitly assumed that the changes of density 
and pressure are small compared to the static values, 
furthermore, the particle veloctiy should be much smaller 
than the sound velocity of the medium. 


Assuming the medium is an ideal gas 


where Y is the ratio of the specific heats for the gas. The 


rd 


speed of sound in a gas is C= ae D.1(a) 
Oo 


For free air the compressiblity and rarification of the 


C 
wave occurs adiabatically and accordingly Son atl 
V 
C = 343 m/s . In small pipes at low frequencies 


(Ks4) the compressibilty is closer to isothermal, thus 
decreasing the propagation speed y= 1, C = 290 m/s ; 
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Substituting (D2) into (D3) yields 


Using again the harmonic variation in time of the 


variables and the mean velocity, the force balance is 
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and equation (D4) becomes 


otha sh oe P D.6 
oC 


taking divergence of (D4) and substituting (D5) yields 
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Applying Galerkin's criterion to (D10) gives 
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[s] {P,} o Neue Leas + ie) bay = 0 Oweiky) 


The boundary condition whe 


aX 16 will be evaluated next. 


For a rigid terminus dP/dX=0 at X=L. The evaluation of 
the boundary condition Nise at X=0 (at the open end of the 


pipe) begins with 
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[D] = {N} LNG D246 


If a linear interpolation function is used 
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The complete equation to be solved is 
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This gives the pressure distribution along the pipe. To 
find the impedence of the pipe at the open end the pressure 
is divided by the velocity at the open end (the velocity is 
that used to calculate {Q}). 

To find the pressure distribution in a pipe of varying 
area a Similar analysis is done by using the wave equation 
fora. horn, which is 
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The only difference in the resulting equations is the 
manner [S] and [P] are calculated and the input term is 


multiplied by the area at the open end. 
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If a linear interpolation function is used 


3A, + A A. +A 
one 1 2 1 2 peo2 
2 2 


ah “grniteeipa inh shure a eee 


1: wee tenet tlie Ko hla wl TAY ba 
Aina vegies ote oe ‘oe ie ef 


: = 4 
- 
cs ee Aa 
Y -*, 


: : Cy Die 1 Se 
vith, gel = on! petal eel Ogos leg: 


7 amie ees’. 
Poa. wae * 


7 aay 


Onde Sept EB LY 
Tee ak Yel ae CK lie chee 


m= Be ae Leek 
Pere ie Sewn TRAE 
ea 


